*** Density ********************************************************************

// Load the data.
use "$pathD/all_tract.dta", clear

// Find the averages.
sum(PM_2p5_indoor) if med_smokepm_pred == 0
local av_ind_nosmoke = r(mean)
local av_ind_nosmoke_1 = r(mean) + .5
local av_ind_nosmoke_num = string(`av_ind_nosmoke', "%4.1f")
sum(PM_2p5_outdoor) if med_smokepm_pred == 0
local av_outd_nosmoke = r(mean)
local av_outd_nosmoke_1 = r(mean) + .5
local av_outd_nosmoke_num = string(`av_outd_nosmoke', "%4.1f")
sum(PM_2p5_indoor) if med_smokepm_pred == 1
local av_ind_lowsmoke = r(mean)
local av_ind_lowsmoke_1 = r(mean) + .5
local av_ind_lowsmoke_num = string(`av_ind_lowsmoke', "%4.1f")
sum(PM_2p5_outdoor) if med_smokepm_pred == 1
local av_outd_lowsmoke = r(mean)
local av_outd_lowsmoke_1 = r(mean) + .5
local av_outd_lowsmoke_num = string(`av_outd_lowsmoke', "%4.1f")
sum(PM_2p5_indoor) if med_smokepm_pred == 2
local av_ind_highsmoke = r(mean)
local av_ind_highsmoke_1 = r(mean) + .5
local av_ind_highsmoke_num = string(`av_ind_highsmoke', "%4.1f")
sum(PM_2p5_outdoor) if med_smokepm_pred == 2
local av_outd_highsmoke = r(mean)
local av_outd_highsmoke_1 = r(mean) + .5
local av_outd_highsmoke_num = string(`av_outd_highsmoke', "%4.1f")

// Generate the needed lines.
sort med_smokepm_pred tract_id date
by med_smokepm_pred: gen x_ind = -1 if _n <= 2
by med_smokepm_pred: gen y_ind = 0 if _n == 1
by med_smokepm_pred: replace y_ind = -2 if _n == 2
replace x_ind = `av_ind_nosmoke' if med_smokepm_pred == 0 & !missing(x_ind)
replace x_ind = `av_ind_lowsmoke' if med_smokepm_pred == 1 & !missing(x_ind)
replace x_ind = `av_ind_highsmoke' if med_smokepm_pred == 2 & !missing(x_ind)
replace y_ind = .25 if y_ind == -2 & med_smokepm_pred == 0
replace y_ind = .2 if y_ind == -2 & med_smokepm_pred == 1
replace y_ind = .075 if y_ind == -2 & med_smokepm_pred == 2
by med_smokepm_pred: gen x_outd = -1 if _n <= 2
by med_smokepm_pred: gen y_outd = 0 if _n == 1
by med_smokepm_pred: replace y_outd = -2 if _n == 2
replace x_outd = `av_outd_nosmoke' if med_smokepm_pred == 0 & !missing(x_outd)
replace x_outd = `av_outd_lowsmoke' if med_smokepm_pred == 1 & !missing(x_outd)
replace x_outd = `av_outd_highsmoke' if med_smokepm_pred == 2 & !missing(x_outd)
replace y_outd = .23 if y_outd == -2 & med_smokepm_pred == 0
replace y_outd = .184 if y_outd == -2 & med_smokepm_pred == 1
replace y_outd = .075 if y_outd == -2 & med_smokepm_pred == 2
gen x = 110 if _n == 1
gen y = 0 if _n == 1
drop if PM_2p5_indoor > 110
drop if PM_2p5_outdoor > 110

// Generate PA figures.
twoway (kdensity PM_2p5_indoor if med_smokepm_pred == 0, lcolor(gs3) ///
		lpattern(solid) lwidth(medthick)) ///
	(kdensity PM_2p5_outdoor if med_smokepm_pred == 0, lcolor(gs3) ///
		lpattern(dash) lwidth(medthick)) ///
	(line y_ind x_ind if med_smokepm_pred == 0, ///
		lwidth(vthin) lcolor(gs9) lpattern(solid)) ///
	(line y_outd x_outd if med_smokepm_pred == 0, ///
		lwidth(vthin) lcolor(gs9) lpattern(solid)) ///
	(scatter y x) ///
	, ///
	yscale(range(0 .25) lcolor(gs3)) ///
	ylabel(0(.05).25, nogrid valuelabel labsize(vsmall) labcolor(gs3) ///
		tlcolor(gs3)) ///
	ytitle("") ///
	xscale(range(0 100) lcolor(gs3)) ///
	xlabel(0(10)40, nogrid valuelabel labsize(vsmall) labcolor(gs3) ///
		tlcolor(gs3)) ///
	xtitle("") ///
	legend(off) ///
	graphregion(margin(10 0 45 0) fcolor(white) ///
		lcolor(white) ifcolor(white) ilcolor(white)) ///
	text(.125 -8 "Density", orient(vertical) ///
		color(gs3) size(small) j(center) placement(c)) ///
	text(.21 -.5 "Indoor", ///
		color(gs3) size(vsmall) j(left) placement(r)) ///
	text(.030 12 "Outdoor", ///
		color(gs3) size(vsmall) j(left) placement(r)) ///
	text(.25 `av_ind_nosmoke_1' ///
		"Av. indoor = `av_ind_nosmoke_num' {&mu}g/m{superscript:3}", ///
		color(gs3) size(vsmall) j(left) placement(r)) ///
	text(.23 `av_outd_nosmoke_1' ///
		"Av. outdoor = `av_outd_nosmoke_num' {&mu}g/m{superscript:3}", ///
		color(gs3) size(vsmall) j(left) placement(r)) ///
	text(-.05 40 "PM{subscript:2.5} ({&mu}g/m{superscript:3})", ///
		color(gs3) size(small) j(right) placement(l))
graph export "$pathR/Figures/dens_nosmoke.pdf", replace

twoway (kdensity PM_2p5_indoor if med_smokepm_pred == 1, lcolor(gs3) ///
		lpattern(solid) lwidth(medthick)) ///
	(kdensity PM_2p5_outdoor if med_smokepm_pred == 1, lcolor(gs3) ///
		lpattern(dash) lwidth(medthick)) ///
	(line y_ind x_ind if med_smokepm_pred == 1, ///
		lwidth(vthin) lcolor(gs9) lpattern(solid)) ///
	(line y_outd x_outd if med_smokepm_pred == 1, ///
		lwidth(vthin) lcolor(gs9) lpattern(solid)) ///
	(scatter y x) ///
	, ///
	yscale(range(0 .20) lcolor(gs3)) ///
	ylabel(0(.05).20, nogrid valuelabel labsize(vsmall) labcolor(gs3) ///
		tlcolor(gs3)) ///
	ytitle("") ///
	xscale(range(0 100) lcolor(gs3)) ///
	xlabel(0(10)40, nogrid valuelabel labsize(vsmall) labcolor(gs3) ///
		tlcolor(gs3)) ///
	xtitle("") ///
	legend(off) ///
	legend(order(1 "Indoor" 2 "Outdoor") ///
		region(lcolor(gs3) lwidth(thin)) color(gs3) size(vsmall) rows(2) ///
		position(3) ring(0)) ///
	graphregion(margin(10 0 45 0) fcolor(white) ///
		lcolor(white) ifcolor(white) ilcolor(white)) ///
	text(.10 -8 "Density", orient(vertical) ///
		color(gs3) size(small) j(center) placement(c)) ///
	text(.1525 0 "Indoor", ///
		color(gs3) size(vsmall) j(left) placement(r)) ///
	text(.0225 22 "Outdoor", ///
		color(gs3) size(vsmall) j(left) placement(r)) ///
	text(.2 `av_ind_lowsmoke_1' ///
		"Av. indoor = `av_ind_lowsmoke_num' {&mu}g/m{superscript:3}", ///
		color(gs3) size(vsmall) j(left) placement(r)) ///
	text(.184 `av_outd_lowsmoke_1' ///
		"Av. outdoor = `av_outd_lowsmoke_num' {&mu}g/m{superscript:3}", ///
		color(gs3) size(vsmall) j(left) placement(r)) ///
	text(-.04 40 "PM{subscript:2.5} ({&mu}g/m{superscript:3})", ///
		color(gs3) size(small) j(right) placement(l))
graph export "$pathR/Figures/dens_lowsmoke.pdf", replace

twoway (kdensity PM_2p5_indoor if med_smokepm_pred == 2, lcolor(gs3) ///
		lpattern(solid) lwidth(medthick)) ///
	(kdensity PM_2p5_outdoor if med_smokepm_pred == 2, lcolor(gs3) ///
		lpattern(dash) lwidth(medthick)) ///
	(line y_ind x_ind if med_smokepm_pred == 2, ///
		lwidth(vthin) lcolor(gs9) lpattern(solid)) ///
	(line y_outd x_outd if med_smokepm_pred == 2, ///
		lwidth(vthin) lcolor(gs9) lpattern(solid)) ///
	(scatter y x) ///
	, ///
	yscale(range(0 .075) lcolor(gs3)) ///
	ylabel(0(.025).075, nogrid valuelabel labsize(vsmall) labcolor(gs3) ///
		tlcolor(gs3)) ///
	ytitle("") ///
	xscale(range(0 25) lcolor(gs3)) ///
	xlabel(0(10)100, nogrid valuelabel labsize(vsmall) labcolor(gs3) ///
		tlcolor(gs3)) ///
	xtitle("") ///
	legend(off) ///
	graphregion(margin(10 0 45 0) fcolor(white) ///
		lcolor(white) ifcolor(white) ilcolor(white)) ///
	text(.0385 -8 "Density", orient(vertical) ///
		color(gs3) size(small) j(center) placement(c)) ///
	text(.06 5 "Indoor", ///
		color(gs3) size(vsmall) j(left) placement(r)) ///
	text(.0225 30 "Outdoor", ///
		color(gs3) size(vsmall) j(left) placement(r)) ///
	text(.075 `av_ind_highsmoke_1' ///
		"Av. indoor = `av_ind_highsmoke_num' {&mu}g/m{superscript:3}", ///
		color(gs3) size(vsmall) j(left) placement(r)) ///
	text(.075 `av_outd_highsmoke_1' ///
		"Av. outdoor = `av_outd_highsmoke_num' {&mu}g/m{superscript:3}", ///
		color(gs3) size(vsmall) j(left) placement(r)) ///
	text(-0.015 100 "PM{subscript:2.5} ({&mu}g/m{superscript:3})", ///
		color(gs3) size(small) j(right) placement(l))
graph export "$pathR/Figures/dens_highsmoke.pdf", replace

*** Deciles of wildfires smoke *************************************************

// Load the data.
use "$pathD/all_tract.dta", clear

// Tab the cut-offs.
bysort c_smokepm_pred: egen max_smokepm_pred = max(smokepm_pred)
tab max_smokepm_pred

forvalues s = 0/10 {
	forvalues i = 0/1 {
		foreach inout in in out {
			qui ci means PM_2p5_`inout'door ///
				if c_smokepm_pred == `s' & above_med_income == `i'
			local `inout'_m_s`s'i`i' = r(mean)
			local `inout'_lb_s`s'i`i' = r(lb)
			local `inout'_ub_s`s'i`i' = r(ub)
			}
		}
	}
qui matrix A = 	0, 0, `in_m_s0i0', `in_lb_s0i0', `in_ub_s0i0', ///
						`out_m_s0i0', `out_lb_s0i0', `out_ub_s0i0' \ ///
				0, 1, `in_m_s0i1', `in_lb_s0i1', `in_ub_s0i1', ///
						`out_m_s0i1', `out_lb_s0i1', `out_ub_s0i1' \ ///
				1, 0, `in_m_s1i0', `in_lb_s1i0', `in_ub_s1i0', ///
						`out_m_s1i0', `out_lb_s1i0', `out_ub_s1i0' \ ///
				1, 1, `in_m_s1i1', `in_lb_s1i1', `in_ub_s1i1', ///
						`out_m_s1i1', `out_lb_s1i1', `out_ub_s1i1' \ ///
				2, 0, `in_m_s2i0', `in_lb_s2i0', `in_ub_s2i0', ///
						`out_m_s2i0', `out_lb_s2i0', `out_ub_s2i0' \ ///
				2, 1, `in_m_s2i1', `in_lb_s2i1', `in_ub_s2i1', ///
						`out_m_s2i1', `out_lb_s2i1', `out_ub_s2i1' \ ///
				3, 0, `in_m_s3i0', `in_lb_s3i0', `in_ub_s3i0', ///
						`out_m_s3i0', `out_lb_s3i0', `out_ub_s3i0' \ ///
				3, 1, `in_m_s3i1', `in_lb_s3i1', `in_ub_s3i1', ///
						`out_m_s3i1', `out_lb_s3i1', `out_ub_s3i1' \ ///
				4, 0, `in_m_s4i0', `in_lb_s4i0', `in_ub_s4i0', ///
						`out_m_s4i0', `out_lb_s4i0', `out_ub_s4i0' \ ///
				4, 1, `in_m_s4i1', `in_lb_s4i1', `in_ub_s4i1', ///
						`out_m_s4i1', `out_lb_s4i1', `out_ub_s4i1' \ ///
				5, 0, `in_m_s5i0', `in_lb_s5i0', `in_ub_s5i0', ///
						`out_m_s5i0', `out_lb_s5i0', `out_ub_s5i0' \ ///
				5, 1, `in_m_s5i1', `in_lb_s5i1', `in_ub_s5i1', ///
						`out_m_s5i1', `out_lb_s5i1', `out_ub_s5i1' \ ///
				6, 0, `in_m_s6i0', `in_lb_s6i0', `in_ub_s6i0', ///
						`out_m_s6i0', `out_lb_s6i0', `out_ub_s6i0' \ ///
				6, 1, `in_m_s6i1', `in_lb_s6i1', `in_ub_s6i1', ///
						`out_m_s6i1', `out_lb_s6i1', `out_ub_s6i1' \ ///
				7, 0, `in_m_s7i0', `in_lb_s7i0', `in_ub_s7i0', ///
						`out_m_s7i0', `out_lb_s7i0', `out_ub_s7i0' \ ///
				7, 1, `in_m_s7i1', `in_lb_s7i1', `in_ub_s7i1', ///
						`out_m_s7i1', `out_lb_s7i1', `out_ub_s7i1' \ ///
				8, 0, `in_m_s8i0', `in_lb_s8i0', `in_ub_s8i0', ///
						`out_m_s8i0', `out_lb_s8i0', `out_ub_s8i0' \ ///
				8, 1, `in_m_s8i1', `in_lb_s8i1', `in_ub_s8i1', ///
						`out_m_s8i1', `out_lb_s8i1', `out_ub_s8i1' \ ///
				9, 0, `in_m_s9i0', `in_lb_s9i0', `in_ub_s9i0', ///
						`out_m_s9i0', `out_lb_s9i0', `out_ub_s9i0' \ ///
				9, 1, `in_m_s9i1', `in_lb_s9i1', `in_ub_s9i1', ///
						`out_m_s9i1', `out_lb_s9i1', `out_ub_s9i1' \ ///
				10, 0, `in_m_s10i0', `in_lb_s10i0', `in_ub_s10i0', ///
						`out_m_s10i0', `out_lb_s10i0', `out_ub_s10i0' \ ///
				10, 1, `in_m_s10i1', `in_lb_s10i1', `in_ub_s10i1', ///
						`out_m_s10i1', `out_lb_s10i1', `out_ub_s10i1'
clear
svmat A
rename A1 c_smokepm_pred
rename A2 above_med_income
rename A3 mean_indoor
rename A4 lower_indoor
rename A5 upper_indoor
rename A6 mean_outdoor
rename A7 lower_outdoor
rename A8 upper_outdoor
replace c_smokepm_pred = c_smokepm_pred + .10 if above_med_income
replace c_smokepm_pred = c_smokepm_pred - .10 if !above_med_income

// Generate PA figure.
twoway (scatter mean_indoor c_smokepm_pred if above_med_income == 0, ///
		mcolor(cranberry%100) msize(small) sort) ///
	(scatter mean_indoor c_smokepm_pred if above_med_income == 1, ///
		mcolor(dknavy%100) msize(small) sort) ///
	(scatter mean_outdoor c_smokepm_pred if above_med_income == 0, yaxis(2) ///
		mcolor(cranberry%100) msize(medlarge) msymbol(X) sort) ///
	(scatter mean_outdoor c_smokepm_pred if above_med_income == 1, yaxis(2) ///
		mcolor(dknavy%100) msize(medlarge) msymbol(X) sort) ///
	(rspike upper_indoor lower_indoor c_smokepm_pred ///
		if above_med_income == 0, sort lcolor(cranberry%100)) ///
	(rspike upper_indoor lower_indoor c_smokepm_pred if above_med_income == 1, ///
		sort lcolor(dknavy%100)) ///
	(rspike upper_outdoor lower_outdoor c_smokepm_pred ///
		if above_med_income == 0, yaxis(2) sort lcolor(cranberry%100)) ///
	(rspike upper_outdoor lower_outdoor c_smokepm_pred ///
		if above_med_income == 1, yaxis(2) sort lcolor(dknavy%100)) ///
	, ///
	yscale(range(-40 50) lcolor(gs3)) ///
	yscale(range(0 450) lcolor(gs3) axis(2)) ///
	ylabel(0(10)50, angle(horizontal) labsize(vsmall) labcolor(gs3) ///
		tlcolor(gs3) glcolor(white)) ///
	ylabel(0(50)200, angle(horizontal) labsize(vsmall) labcolor(gs3) ///
		tlcolor(gs3) glcolor(white) axis(2)) ///
	ytitle("") ytitle("", axis(2)) ///
	xscale(range(-.5 10.5) lcolor(gs3)) ///
	xlabel(none) ///
	xticks(0(1)10) ///
	xtitle("") ///
	yline(-10, lcolor(gs9) lwidth(thin) lpattern(dash)) ///
	yline(-20, lcolor(gs9) lwidth(thin) lpattern(dash)) ///
	yline(-30, lcolor(gs9) lwidth(thin) lpattern(dash)) ///
	yline(-40, lcolor(gs3) lwidth(thin) lpattern(dash)) ///
	yline(0, lcolor(gs3) lwidth(thin) lpattern(solid)) ///
	yline(10, lcolor(gs9) lwidth(thin) lpattern(dash)) ///
	yline(20, lcolor(gs9) lwidth(thin) lpattern(dash)) ///
	yline(30, lcolor(gs9) lwidth(thin) lpattern(dash)) ///
	yline(40, lcolor(gs9) lwidth(thin) lpattern(dash)) ///
	yline(50, lcolor(gs9) lwidth(thin) lpattern(dash)) ///
	yline(52.3, lcolor(gs3) lwidth(thin) lpattern(solid)) ///
	legend(order(1 "Low income indoor (left axis)" ///
		2 "High income indoor (left axis)" ///
		5 "Low income outdoor (right axis)" ///
		6 "High income outdoor (right axis)") ///
		region(lcolor(gs9) lwidth(thin)) color(gs3) size(vsmall) rows(4) ///
		position(11) ring(0)) ///
	graphregion(margin(50 5 7 1) fcolor(white) ///
		lcolor(white) ifcolor(white) ilcolor(white)) ///
	text(-48.5 0 "No" "smoke", ///
		color(gs3) size(vsmall) j(c) placement(c)) ///
	text(-48.5 1 "< 1.2", ///
		color(gs3) size(vsmall) j(c) placement(c)) ///
	text(-48.5 2 "1.2" "-" "2.3", ///
		color(gs3) size(vsmall) j(c) placement(c)) ///
	text(-48.5 3 "2.3" "-" "3.5", ///
		color(gs3) size(vsmall) j(c) placement(c)) ///
	text(-48.5 4 "3.5" "-" "5.4", ///
		color(gs3) size(vsmall) j(c) placement(c)) ///
	text(-48.5 5 "5.4" "-" "7.6", ///
		color(gs3) size(vsmall) j(c) placement(c)) ///
	text(-48.5 6 "7.6" "-" "10.1", ///
		color(gs3) size(vsmall) j(c) placement(c)) ///
	text(-48.5 7 "10.1" "-" "13.8", ///
		color(gs3) size(vsmall) j(c) placement(c)) ///
	text(-48.5 8 "13.8" "-" "20.6", ///
		color(gs3) size(vsmall) j(c) placement(c)) ///
	text(-48.5 9 "20.6" "-" "41.9", ///
		color(gs3) size(vsmall) j(c) placement(c)) ///
	text(-48.5 10 "> 41.9", ///
		color(gs3) size(vsmall) j(c) placement(c)) ///
	text(25 -1.9 "Indoor PM{subscript:2.5} ({&mu}g/m{superscript:3})", ///
		color(gs3) size(small) j(left) placement(r) orient(vertical)) ///
	text(100 12 ///
		"Outdoor PM{subscript:2.5} ({&mu}g/m{superscript:3})", color(gs3) ///
			size(small) j(right) placement(l) orient(vertical) yaxis(2)) ///
	text(-66.9 10.5 ///
		"Deciles of wildfire smoke PM{subscript:2.5} ({&mu}g/m{superscript:3})", ///
		color(gs3) size(small) j(right) placement(l))
graph export "$pathR/Figures/wildfire_indoor_outdoor_PM.pdf", replace

clear
